## Difference-in-Differences Analysis (Results are estimated by lfe package version 2.8-7)
## These results are reported in Table 2
rm(list=ls())
options(scipen=100,digits=10)
library(raster)
library(plyr)
library(dplyr)
library(lfe)
library(ggplot2)
library(stargazer)
load("data_final.RData")

########################
# Generate sample without post-shutdown period
data4=subset(data3, T<99)
########################
# Generate the short sample
data5=subset(data4, T>40)

##############################################################################################################################
########################
## Baseline Model
data_base1=data4
data_base2=data5
## Long sample
all_AOD_long=felm(AOD~ shutdown  + ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs.  +factor(Year)+factor(weekday) | factor(weeks)*factor(Plant_Code) + factor(T),
                  data = data_base1)
summary(all_AOD_long)
all_SO2_long=felm(SO2_MASS..tons.~ shutdown  + ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs. +factor(Year)+factor(weekday) | factor(weeks)*factor(Plant_Code) + factor(T),
                  data = data_base1)
summary(all_SO2_long)
all_NOX_long=felm(NOX_MASS..tons.~ shutdown  + ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs. +factor(Year)+factor(weekday) | factor(weeks)*factor(Plant_Code) + factor(T),
                  data = data_base1)
summary(all_NOX_long)

## Short sample
all_AOD_short=felm(AOD~ shutdown + ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs.  +factor(Year)+factor(weekday) | factor(weeks)*factor(Plant_Code) + factor(T),
                   data = data_base2)
summary(all_AOD_short)
all_SO2_short=felm(SO2_MASS..tons.~ shutdown + ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs. +factor(Year)+factor(weekday) | factor(weeks)*factor(Plant_Code) + factor(T),
                   data = data_base2)
summary(all_SO2_short)
all_NOX_short=felm(NOX_MASS..tons.~ shutdown + ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs. +factor(Year)+factor(weekday) | factor(weeks)*factor(Plant_Code) + factor(T),
                   data = data_base2)
summary(all_NOX_short)
## Output results
stargazer(all_SO2_long, all_NOX_long, all_AOD_long,all_SO2_short, all_NOX_short, all_AOD_short)
########################
## Cluster s.d. at plant level
## Long sample
all_AOD_long_C=felm(AOD~ shutdown  + ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs.  +factor(Year)+factor(weekday) | factor(weeks)*factor(Plant_Code) + factor(T)|0|Plant_Code,
                    data = data_base1)
summary(all_AOD_long_C)
all_SO2_long_C=felm(SO2_MASS..tons.~ shutdown  + ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs. +factor(Year)+factor(weekday) | factor(weeks)*factor(Plant_Code) + factor(T)|0|Plant_Code,
                    data = data_base1)
summary(all_SO2_long_C)
all_NOX_long_C=felm(NOX_MASS..tons.~ shutdown  + ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs. +factor(Year)+factor(weekday) | factor(weeks)*factor(Plant_Code) + factor(T)|0|Plant_Code,
                    data = data_base1)
summary(all_NOX_long_C)
## Short sample
all_AOD_short_C=felm(AOD~ shutdown + ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs.  +factor(Year)+factor(weekday) | factor(weeks)*factor(Plant_Code) + factor(T)|0|Plant_Code,
                     data = data_base2)
summary(all_AOD_short_C)
all_SO2_short_C=felm(SO2_MASS..tons.~ shutdown + ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs. +factor(Year)+factor(weekday) | factor(weeks)*factor(Plant_Code) + factor(T)|0|Plant_Code,
                     data = data_base2)
summary(all_SO2_short_C)
all_NOX_short_C=felm(NOX_MASS..tons.~ shutdown + ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs. +factor(Year)+factor(weekday) | factor(weeks)*factor(Plant_Code) + factor(T)|0|Plant_Code,
                     data = data_base2)
summary(all_NOX_short_C)
## Output results
stargazer(all_SO2_long_C, all_NOX_long_C,all_AOD_long_C,all_SO2_short_C, all_NOX_short_C, all_AOD_short_C)
